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Chapter  1 

Introduction  and  Overview 


This  final  report  contains  a  summary  of  the  activities  supported  under  the  Air  Force  AFOSR 
AASERT  grant  F49620-97-1-0356.  This  grant  supported  two  graduate  students,  Lisa  Stanley  and 
Kevin  Hulsing,  to  work  with  Professor  John  A.  Burns  (currently  supported  under  AFOSR  PRET 
Grant  F49620-96- 1-0329)  on  computational  methods  for  optimal  design  and  control  of  aerospace 
systems.  In  addition,  the  grant  provided  partial  support  for  Major  Dawn  Stewart  who  completed 
her  Ph.D  at  Virginia  Tech.  The  research  is  based  on  our  current  work  on  mathematical  and  com¬ 
putational  methods  for  sensitivity  analysis  and  on  new  controller  reduction  methods  for  systems 
governed  by  partial  differential  equations.  This  work  is  fully  described  in  the  proposal  “Sensitivity 
and  Adjoint  Methods  for  Design  of  Aerospace  Systems”  (this  proposal  resulted  in  the  AFOSR  PRET 
grant).  The  focus  of  this  project  is  fundamental  research  in  sensitivity  based  methods  for  optimal 
design  and  in  computational  methods  for  practical  distributed  parameter  control.  A  major  goal  of 
this  effort  is  to  insure  the  transition  of  this  work  into  Air  Force  and  industrial  applications.  The 
research  conducted  by  the  AASERT  students  is  motivated  by  the  need  to  develop  accurate  solvers 
for  sensitivity  equations  and  to  develop  computational  tools  for  controller  reduction  in  2D  and  3D 
physics  based  models.  The  principal  investigator  is  Dr.  John  A.  Burns,  Director  of  the  Center  for 
Optimal  Design  and  Control  and  Hatcher  Professor  of  Mathematics  at  Virginia  Tech. 

The  AA5J5i2nstudents  worked  within  the  PRET  Center  at  Virginia  Tech.  Research  under  the 
PRET  Grant  is  conducted  at  six  institutions:  Virginia  Tech,  Cornell,  AeroSoft,  Beam  Technologies, 
Boeing  and  Lockheed  Martin.  Frequent  meetings  of  the  team,  exchange  of  visits,  sharing  of  software, 
exchange  of  graduate  students  and  postdoctoral  researchers,  industry- Air  Force  laboratory-university 
workshops  and  communication  of  results  are  coordinated  by  the  Center. 

The  PRET  Center  offers  an  unparalleled  potential  for  strengthening  the  educational  and  scien¬ 
tific  infrastructure  by  training  students  and  post-doctoral  researchers  in  an  interdisciplinary  team 
approach  to  scientific  and  engineering  research.  The  Center  provides  unique  opportunities  for  the¬ 
oretical,  computational,  and  experimental  research.  Through  the  interactions  with  Air  Force  labo¬ 
ratories  and  industrial  partners,  students  are  exposed  to  real  problems.  The  combined  theoretical, 
computational,  and  experimental  approach  provides  a  meaningful  interdisciplinary  research  experi¬ 
ence. 

The  AASERT  students  played  an  integral  role  in  the  transition  of  research  into  industry  and  Air 
Force  laboratories. 
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Chapter  2 

Objectives 


This  ASSERT  grant  was  designed  to  support  the  PRET  program  through  the  support  of  graduate 
students  and  their  research.  In  addition  to  conducting  basic  research,  these  students  played  an 
integral  role  in  the  development  of  computational  tools  and  in  the  transitioning  of  these  tools  to 
industry  and  Air  Force  laboratories. 

2.1  Research  Objectives 

Although  the  PRET  proposal  contains  several  research  objectives  and  the  AASERT  students  are 
involved  in  several  aspects  of  these  projects,  the  AASERT  students  focused  on  two  specific  tasks 
with  the  following  research  objectives. 

•  The  development  of  fast  and  accurate  computational  methods  for  sensitivity  anal¬ 
ysis.  Sensitivity  analysis  plays  a  central  role  in  optimal  design  as  well  as  in  the  development 
of  fast  simulation  algorithms.  The  goal  of  this  effort  was  to  improve  the  accuracy  and  effi¬ 
ciency  in  sensitivity  computations  by  using  adaptivity  and  smoothing  methods.  In  addition, 
the  method  of  mappings  was  investigated  in  order  to  improve  accuracy. 

•  The  development  of  computational  tools  specifically  for  control  design  of  dis¬ 
tributed  parameter  systems.  The  goal  here  is  to  develop  practical  computational  tools 
for  sensor/actuator  placement  and  reduced  order  controller  design.  This  goal  was  achieved  by 
using  distributed  parameter  theory,  Chandrasekhar  partial  differential  equations,  and  repre¬ 
sentations  as  a  foundation  for  computational  algorithms. 

2.2  Educational  Objectives,  Interactions  and  Transitions 

The  educational  goal  of  this  ASSERT  was  to  train  interdisciplinary  applied  mathematicians,  engi¬ 
neers  and  scientists  to  work  on  industrial  problems  of  importance  to  the  Air  Force.  This  goal  was 
met  through  interactions  with  PRET  partners  and  the  corresponding  transition  of  research  into 
these  industries.  In  particular,  we  had  the  following  specific  educational  objectives: 

•  To  develop  an  active  program  directed  at  educating  interdisciplinary  scientists  and 
encouraging  these  graduates  to  work  at  Air  Force  facilities,  laboratories  and  in¬ 
dustry.  All  of  the  students  finished  their  Ph.D.  degrees  and  are  working  on  research  important 
to  the  Air  Force. 

•  To  transition  research  and  computational  tools  into  industry  and  Air  Force  labo¬ 
ratories.  We  worked  with  our  PRET  partners  to  test  new  algorithms  on  design  and  control 
applications  including  (but  not  limited  to)  shape  optimization  and  control  design  for  flow 
management  and  heat  transfer. 
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Chapter  3 

Status  of  Effort 


Considerable  progress  was  made  toward  achieving  the  all  of  the  proposed  research  objectives.  The 
students  developed  new  computational  tools  and  began  the  construction  of  a  mathematical  frame¬ 
work  that  will  allow  for  the  analysis  of  these  tools.  This  effort  focused  on  the  difficulties  of  producing 
practical  design  and  control  tools  for  complex  2D  and  3D  problems.  This  report  contains  a  sum¬ 
mary  of  new  findings  and  accomplishments  (Chapter  4)  and  the  interact  ions/ transitions  (Chapter 
7).  However,  the  following  items  are  particularly  noteworthy. 


Status  of  Effort 

During  this  period,  the  two  students  supported  under  this  grant  produced  7  scientific  papers, 
given  more  than  10  invited  presentations  and  developed  preliminary  software  packages  for  design  and 
control.  In  addition,  the  students  made  considerable  progress  on  three  specific  projects:  adaptive 
finite  element  and  projection  methods  for  sensitivity  computations,  mapping  techniques  for  sensitiv¬ 
ity  computations  in  elliptic  systems,  and  numerical  algorithms  for  computation  of  functional  gains 
for  feedback  control  of  thermal  processes. 

Ms.  Lisa  Stanley,  working  with  Major  Dawn  Stewart  (USAF),  developed  a  new  algorithm  for 
sensitivity  calculations  that  is  based  on  smoothing  projections  of  spatial  derivatives.  They  considered 
both  local  and  global  projections  and  tested  the  numerical  methods  on  several  boundary  value 
problems.  This  work  concentrated  on  projection  error  techniques  for  both  state  and  sensitivity 
equations. 

Major  Stewart  extended  these  methods  to  the  Navier-Stokes  equations  and  enhanced  the  per¬ 
formance  of  the  basic  algorithm  by  employing  adaptive  mesh  generators  to  improve  accuracy.  Ms. 
Stanley  investigated  the  method  of  mappings  and  used  the  continuous  sensitivity  method  to  produce 
accurate  computational  schemes.  She  showed  that,  if  one  first  derives  the  sensitivity  equation  in 
the  physical  domain  and  then  transforms  both  state  and  sensitivity  equation  to  the  computational 
domain,  then  mesh  gradients  are  not  needed  and  the  accuracy  of  computed  sensitivities  can  be  con¬ 
trolled.  This  work  also  provides  a  theoretical  framework  for  the  rigorous  analysis  of  the  resulting 
algorithms. 

Mr.  Kevin  Hulsing  developed  a  fast  new  approach  to  the  problem  of  computing  functional  gains 
that  define  optimal  feedback  controllers  for  spatially  distributed  systems.  He  derived  Chandrasekhar 
partial  differential  equations  for  the  gains  and  constructed  numerical  methods  for  solving  these 
systems.  He  applied  this  method  to  ID  and  2D  thermal  control  problems.  This  approach  is  based 
on  direct  solutions  of  Riccati  and  Chandrasekhar  partial  differential  equations  and  allows  the  use  of 
adaptive  finite  element  methods  (both  p-refinement  and  h-refinement).  In  addition,  he  made  use  of 
“Nitsche”  finite  element  techniques  to  enhance  accuracy  near  the  boundaries. 
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Chapter  4 

Accomplishments  and  New 
Findings 


This  chapter  contains  a  detailed  summary  of  some  of  the  research  accomplishments  and  new  findings, 
a  brief  description  of  research  highlights  and  their  significance  to  the  field.  In  addition,  we  indicate 
the  relevance  of  the  research  to  potential  Air  Force  applications  and  industrial  technology  challenges. 
In  §  4.1  we  present  the  highlights  of  Ms.  Stanley’s  work  on  the  use  of  transformations  in  sensitivity 
computations  for  elliptic  systems.  In  §  4.2  we  discuss  Mr.  Hulsing’s  work  on  numerical  methods  for 
computing  functionsl  gains  for  distributed  processes.  The  following  accomplishments  are  particularly 
important. 

Accomplishment  1:  Developed  a  Method  of  Mappings  Method  for  Sensitivity  Compu¬ 
tations  in  Elliptic  Boundary  Value  Problems 

•  New  Findings:  This  algorithm  is  based  on  first  deriving  the  sensitivity  equations  in  the 
physical  domain  and  then  mapping  both  state  and  sensitivity  equations  to  the  computational 
domain.  This  approach  allows  for  greater  control  of  errors  and  eliminates  the  need  to  com¬ 
pute  mesh  gradients.  In  particular,  these  results  show  that  errors  in  computing  geometry 
sensitivities  are  important,  but  can  be  controlled  by  improving  the  boundary  approximations 
alone. 

•  Significance  and  Potential  Applications:  Computing  sensitivities  in  aerospace  systems 
is  fundamental  to  any  design  process.  For  large  systems  of  differential  equations,  speed  and 
accuracy  of  computation  are  essential.  This  work  provides  insight  into  the  types  of  errors  that 
must  be  controlled  in  order  to  achieve  accuracy.  In  addition,  the  ability  to  compute  sensitivity 
flow  fields  in  a  fraction  of  the  time  it  takes  to  solve  the  basic  equations,  reduces  design  cycle 
times  and  allows  the  designer  to  evaluate  the  model  as  well  as  the  design.  These  types  of 
problems  occur  in  several  Air  Force  applications.  The  design  and  analysis  of  the  COIL  laser 
is  just  one  of  many  such  applications  we  are  considering. 

• 

Accomplishment  2:  Developed  a  Chandrasekhar  PDE  Approach  for  Computing  Func¬ 
tional  Gains  for  Feedback  Control  of  Thermal  Processes 

•  New  Findings:  This  method  is  based  on  deriving  a  system  of  Chandrasekhar  type  partial 
differential  equations  that  can  be  solved  for  the  functional  gains  that  define  optimal  feedback 
laws.  The  motivation  is  to  develop  new  computational  algorithms  that  can  be  applied  to 
2D  and  3D  distributed  parameter  systems.  High  order  accurate  schemes  are  used  to  solve 
the  resulting  system  of  partial  differential  equations.  This  method  was  applied  to  a  2D  heat 
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transfer  problem.  Even  when  all  other  standard  methods  failed,  the  Chandrasekhar  method 
converged  to  the  gains.  Therefore,  this  work  established  the  feasibility  of  the  approach  for 
complex  physisc  based  models. 

•  Significance  and  Potential  Applications:  In  order  to  use  functional  gains  for  sensor/actuator 
placement  and  reduced  order  controller  design,  one  first  needs  accurate  approximations  to 
these  kernel  functions.  This  new  method  not  only  improves  accuracy,  but  also  can  provide 
this  accuracy  on  coarse  grids.  In  addition,  the  algorithm  is  highly  parallelizable  and  allows  for 
adaptivity.  Thus,  this  type  of  method  is  more  suitable  for  2D  and  3D  control  problems  that 
occur  in  flow  control  and  control  of  manufacturing  processes. 
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4.1  Hybrid  Transformation  Methods  for  Sensitivity  Compu¬ 
tations 

Sensitivity  analysis  is  an  important  tool  in  engineering  with  applications  ranging  from  optimal  design 
to  computational  fluid  dynamics.  There  are  many  algorithms  available  for  computing  sensitivities. 
These  range  from  automatic  differentiation  techniques  to  variational  methods.  In  this  section  we 
describe  two  variational  methods  for  computing  sensitivities.  We  use  a  ID  model  problem  as  a 
platform  for  describing  and  comparing  these  methods.  In  particular,  the  regularity  of  the  sensitivity 
equations  obtained  from  each  method  is  addressed.  Numerical  approximations  to  each  sensitivity 
are  calculated  using  finite  element  schemes,  and  numerical  comparisons  are  presented. 

4.1.1  Introduction 

Accurate  sensitivity  calculations  play  an  important  role  in  the  analysis  and  optimization  of  engi¬ 
neering  systems.  Sensitivities  can  be  used  to  compute  gradients  in  optimization-based  design.  In 
addition,  they  have  been  used  to  construct  fast  solvers  for  computational  fluid  dynamics.  In  this  sec¬ 
tion,  we  focus  on  variational  methods  for  computing  state  sensitivities.  These  schemes  make  use  of 
sensitivity  equation  methods.  However,  there  are  a  variety  of  ways  to  implement  sensitivity  equation 
methods,  and  these  variations  yield  algorithms  with  different  convergence  properties.  We  consider 
two  specific  methods.  The  first  is  based  on  transforming  the  state  equation  to  a  fixed  computational 
domain  and  then  deriving  its  sensitivity  equation.  Once  the  state  and  sensitivity  systems  are  solved, 
the  solutions  are  mapped  back  to  the  physical  domain.  The  second  approach  transforms  both  state 
and  sensitivity  equations,  solves  the  transformed  equations  and  maps  these  solutions  back  to  the 
physical  domain.  There  axe  benefits  and  drawbacks  to  each  method.  Indeed,  it  is  not  always  obvious 
which  scheme  is  best  for  a  given  problem.  Many  questions  need  to  be  addressed  before  a  complete 
theory  can  be  developed.  In  the  discussion  below,  we  will  describe  a  simple  ID  example  in  order  to 
illustrate  the  methods  and  to  motivate  some  of  these  questions.  We  begin  with  a  description  of  the 
model  problem.  Each  of  the  variational  methods  will  be  described  in  the  context  of  the  example, 
and  numerical  results  will  be  presented.  Finally,  we  will  compare  the  performance  of  these  methods 
and  make  some  concluding  remarks  concerning  future  work. 


4.1.2  Notation 


Before  describing  the  variational  methods,  we  define  the  transformations  used  to  move  between  the 
physical  and  the  computational  domains.  Some  notation  regarding  Sobolev  spaces  is  also  introduced. 
The  “physical  domain”  for  the  model  problem  is  the  interval  (0,g),  where  q  is  a  parameter  with 
q  €  (1,2).  The  computational  domain  is  the  unit  interval  (0,1).  For  a  >  0,  let  Qa  =  (0,a),  and 
for  each  fixed  q  €  (1,2)  define  the  transformation  T  :  Qq  -4  fli  by  T(x,q)  =  |  =  £.  Note  that  the 
function  M:  fli  -4  Clq  defined  by  M(f)  =  £q  =  x  is  the  inverse  of  T  and  is  commonly  referred  to  as 
the  “mesh  map” . 

Let  Hm(Sl  i)  denote  the  usual  Sobolev  space  of  “functions”  whose  partial  derivatives,  up  to  order 
m,  axe  also  square  integrable.  Let  L 2  =  L2(Q i)  with  inner  product  defined  by 

(u,v)=  f  )v(Od£ 

Jo 

for  all  €  L2.  For  this  study,  we  need  only  to  consider  V  =  Hq(Q i)  c  iT1  (fii ),  where  V 

consists  of  functions  in  H1^ i)  with  zero  trace.  The  dual  space  of  V  is  given  by  V *  =  #-1( fli). 
Finally,  we  need  the  bilinear  form  o:  V  x  V  -*  E  by 

<4>,v)  =  (4.1.1) 

Jo 


for  all  <)>(■),  v(-)  e  V. 
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4.1.3  A  Model  Problem 


In  this  section,  a  ID  model  problem  is  introduced.  This  example  is  convenient  for  numerical  com¬ 
parisons  since  both  the  state  and  the  sensitivity  equations  have  analytical  solutions.  The  model 
problem  is  structured  so  that  the  physical  domain  is  parameter  dependent. 

Let  1  <  q  <  2,  and  consider  the  state  equation  given  by  the  following  elliptic  boundary  value 
problem 


-w"(x)  =  f(x),  x€(0,q) 

(4.1.2) 

with  boundary  conditions 

tu(0)  =  0,  w(q)  =  0. 

(4.1.3) 

Here  /  :  (0,  +oo)  — >  M  is  the  piecewise  continuous  function  given  by 

f(x)  _/  °>  0<x<l 

n  ’  {  -1,  1  <  a;  <  +oo. 

(4.1.4) 

The  parameter  q  determines  the  length  of  the  interval  over  which  the  state  is  defined.  The  goals  are 
to  solve  (4.1.2)-(4.1.3)  for  the  state,  w(x,q)1  for  a  given  value  of  q  and  to  determine  the  sensitivity 
of  the  state  to  small  changes  in  the  parameter.  The  sensitivity  is  denoted 

s(x,q)  =  dqw(x,q). 

(4.1.5) 

For  this  example,  the  state  and  sensitivity  are  given  by 

.  .  f  0  <  x  <  1 

l  J2 q  *+|(*  1)2>  1  <  x  <  q 

(4.1.6) 

and 

s(x,q)=  ^  2q2  ^  x,  xe(0 ,q), 

(4.1.7) 

respectively.  However,  we  proceed  to  solve  for  both  the  state  and  the  sensitivity  numerically  using 
each  of  the  methods  previously  mentioned.  Comparing  the  numerical  calculations  to  (4.1.6)-(4.1.7) 
allows  one  to  quantitatively  measure  the  usefulness  and  accuracy  of  the  methods.  Note  that  for  a 
fixed  q  >  1,  the  function  x  — >  w(x,q)  belongs  to  H2(0,q)  PI  #o(0,g).  However,  dlw(x,  q)  is  discon¬ 
tinuous,  while  the  sensitivity  function,  x  ->  s(x,q)  is  C°°.  Having  defined  the  relevant  equations 
and  unknowns,  we  will  now  use  (4.1.2)-(4.1.3)  to  describe  two  methods  for  computing  numerical 
approximations  to  s(x,q). 

4.1.4  Methods  for  Computing  the  Sensitivity 

For  many  engineering  applications,  a  typical  approach  to  such  problems  is  to  begin  by  transform¬ 
ing  the  problem  to  a  fixed  “computational  domain” .  This  “computational  domain”  is  often  more 
regular  in  shape  which  simplifies  grid  generation  and  can  improve  the  accuracy  of  numerical  calcu¬ 
lations.  This  technique  is  especially  common  for  problems  in  computational  fluid  dynamics  and  for 
many  problems  involving  moving  boundaries  such  as  shape  optimization.  The  transformations  can 
be  performed  at  the  infinite  dimensional  level  as  long  as  the  mappings  are  isomorphisms  and  are 
sufficiently  smooth  in  comparison  to  the  regularity  of  the  solution  to  the  PDE  in  question.  Both 
methods  considered  here  make  use  of  these  transformations.  It  is  important  to  note  that  in  2D 
and  3D  problems,  transforming  can  be  a  complex  process.  In  particular,  in  addition  to  problems 
of  regularity,  it  is  important  to  consider  accuracy  loss.  For  example,  transformations  that  produce 
coordinate  systems  that  are  orthogonal  at  boundaries  are  preferred.  Otherwise,  accuracy  decreases 
as  orthogonality  declines.  For  the  ID  example  here,  it  is  straightforward.  Some  of  the  difficulties 
that  occur  in  2D  and  3D  problems  are  not  present  in  this  case.  In  particular,  this  step  often  requires 
numerical  approximations  of  the  mapping  and  its  spatial  derivatives.  More  about  this  issue  will  be 
discussed  in  the  conclusion. 
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4.1.5  Two  Approaches  for  Sensitivity  Calculation 

Beginning  with  the  state  equation  defined  on  the  physical  domain,  one  can  implicitly  differentiate 
the  state  (4.1.2)-(4.1.3)  in  order  to  obtain  a  sensitivity  equation.  At  this  stage,  it  is  important 
to  note  that  this  differentiation  is  rather  formal.  In  general,  the  partial  derivatives  d*w{x,  q )  and 
dqw(x ,  q)  need  to  be  continuous  in  order  to  interchange  the  order  of  differentiation.  In  this  case 
dlw(x,q)  is  discontinuous,  but  one  can  verify  by  hand  that  the  sensitivity  dqw(x,q)  satisfies  the 
following  differential  equation 


-s"(x)  =0,  0  <  x  <  q 


(4.1.8) 


with  boundary  conditions 


s(0)  =  0,  s(q)  =  -dxw(q)  ~  -dxw(x) 


(4.1.9) 


x=g 


System  (4.1.8)-(4.1.9)  is  called  the  sensitivity  equation.  The  boundary  conditions  should  be  derived 
with  care  as  the  right  endpoint  of  the  domain  (0,  q)  depends  explicitly  on  the  parameter  g. 

Once  the  sensitivity  equation  has  been  derived  on  the  physical  domain,  the  transformations  in 
Section  4.1.2  are  used  to  define  the  “transformed”  functions.  For  £  E  fii  and  q  E  (1,2),  define 


w(^q)  =w(M(£,q)iq)  =  w(x9q), 
*(f»9)  =  s(Af(f,g),g)  =  s(x9q) 


and 


/(£>«)  =  f(M(Z,q),q)  -  f(x). 
Once  transformed,  the  original  forcing  function  /  becomes 


S(C  \  _  /  °>  0<?<i 

/(e,9)  1-1,  I<^<1 


(4.1.10) 

(4.1.11) 

(4.1.12) 


and  now  depends  explicitly  on  the  parameter  q.  Using  the  above  definitions  and  the  chain  rule,  the 
spatial  derivatives  of  the  original  functions  and  those  of  the  transformed  functions  are  related  by 

dxw(x,q)  =  diw(T(x,q),q)-dxT(x,q) 

=  (4.1.13) 

and 

d%w(x,q)  =  d\w(T{x,q),q)  -[dxT{x,q)]2  + 

diw(T(x,q),q)-dlT(x,q ) 

=  d2(w(€,q)'^-  (4.1.14) 

These  identities  are  used  to  derive  transformed  boundary  value  problems  for  both  the  transformed 
state  and  the  transformed  sensitivity  in  (4.1.10).  The  transformed  state  equation  is  defined  by  the 
differential  equation 


with  boundary  conditions 


-*"(£)  =  «2 /(£,?),  fe(o,i), 
w(0)  =  0,  u>(l)  =  0. 


(4.1.15) 

(4.1.16) 
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Likewise,  the  transformed  sensitivity  equation  is  given  by  the  differential  equation 


-S"(f)=0,  £€(0,1), 


(4.1.17) 


with  boundary  conditions 


5(0)  =  0, 


•W)  =  -(i) 


•  %«>(£) 


l€=i 


(4.1.18) 


Using  (4.1.6)  and  (4.1.7),  the  solutions  to  the  transformed  state  and  sensitivity  equations  are  given 
by 


and 


«>(£.?)  = 


{ 


.k 


ill! 


2 

(9~1): 

2 


£, 

£  + 


o<£<^ 


(4.1.19) 


*(£,  ff)  =  --?2„  -  o  <  £  <  1,  (4.1.20) 

2  q 

respectively. 

Before  describing  the  variational  formulations  and  the  discretization,  we  will  discuss  an  alter¬ 
nate  approach.  The  process  we  outline  in  the  following  paragraphs  varies  from  the  one  previously 
described  in  the  order  in  which  the  sensitivity  equation  is  derived  and  the  transformations  are 
performed. 

A  second  approach  to  the  computation  of  the  sensitivity  is  similar  in  spirit  to  the  Semi- Analytical 
method  (SAM).  This  technique  is  often  used  in  the  engineering  community.  Roughly  speaking,  the 
SAM  begins  by  first  transforming  the  state  equation  to  the  computational  domain.  The  second  step 
is  to  discretize  the  state  equation,  thereby  producing  an  algebraic  system.  This  discrete  equation  is 
then  differentiated  to  obtain  a  discrete  sensitivity  equation  which  is  solved  using  special  techniques. 
An  abstract  version  of  this  method  (A-SAM)  may  be  constructed  by  deriving  a  sensitivity  equation 
after  transforming  but  before  discretizing  the  state  equation.  In  particular,  the  infinite  dimensional 
transformed  state  equation  is  “differentiated”  in  order  to  obtain  an  equation  for  the  sensitivity  of 
the  transformed  state.  We  now  present  the  details  of  this  approach. 

On  the  computational  domain,  fix,  define  the  sensitivity  of  the  transformed  state  by  p(£, q)  — 
9qw(^q).  In  order  to  derive  a  system  for  p(£, g),  the  transformed  state  equation  (4.1.15)-(4.1.16)  is 
“differentiated”  with  respect  to  g.  This  requires  that  the  function  /(£,  q)  be  “differentiated”  with 
respect  to  q  as  well.  As  before,  formal  differentiation  yields  the  boundary  value  problem 


-dfp(£.  q)  =  $(£,  q),  (4.1.21) 

p(0)  =  0,  p{  1)  =  0,  (4.1.22) 

where  g(£,q)  =  2qf(£,q)  +  q2dqf(£,  q)  belongs  to  V* .  In  particular, 

0(£,«)  =  2ff/(£,g)-*i(£),  (4.1.23) 

where  S i(f)  is  the  delta  function  with  mass  at  L.  Since  the  linear  elliptic  problem  (4.1.21)-(4.1.22) 

does  not  have  £f2(fli)  solutions,  the  system  must  be  considered  in  the  weak  sense.  We  use  this 
system  to  derive  a  weak  formulation  in  Section  4.1.6. 

For  this  example,  p(£,  q)  can  be  calculated  directly  from  (4.1.19)  to  obtain 


dqw(£,q) 


~(q  - 1)£> 

~(q  - 1)£  +  (£?  - 1)£, 


0<£<i 

?<£<!• 


(4.1.24) 
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Clearly,  dqw(x ,  q )  =  s(£,  q )  ^  p(£,  q).  In  particular,  we  observe  that  the  sensitivity  of  the  transformed 
state,  p(f,g),  is  less  smooth  than  the  transformed  sensitivity  s(£,  q)  in  (4.1.20).  Furthermore,  the 
relationship  between  the  transformed  sensitivity  and  the  sensitivity  of  the  transformed  state  can  be 
derived  using  (4.1.13)-(4.1.14)  and  the  chain  rule.  Beginning  with  the  definition  of  the  sensitivity  of 
the  transformed  state,  it  follows  that 

?(£>?)  =  dqw(Z,  q) 

=  dqw(M(Z,q),q ) 

=  dqw(M(Z,  q),q)  +  dxw(M{Z,  q),  q)  [dqM(Z,  9)] 

=  s(M(£,  q),q)  +  dxw(x,q)  [dqM(Z,q)} 

=  s(£,  q)  +  [dxT{ x,  9)]  9)  [dqM (f,  9)] 

=  s(Z,q)  +  9)]-1  d(w(Z,  9)  [dqM(Z,  9)]  • 

The  relationship  between  s(x,  9)  and  p(Z,  9)  can  be  obtained  by  using  (4.1.10)  and  the  definition  of 
the  transformation  T.  Direct  computation  yields 

s(z,g)  =p(^,q) -dsw(^,q)  |d?M(^,9)j  ^d„M(^,q).  (4.1.25) 

Observe  the  appearance  of  the  “mesh  derivative” ,  dqM(< f ,  g),  and  the  spatial  derivative  of  u)(£,  q)  in 
this  equation.  In  order  to  calculate  s(x,q)  using  this  (A-SAM)  approach,  one  needs  to  compute  not 
only  the  sensitivity  of  the  transformed  state  but  also  the  spatial  derivative  of  the  transformed  state 
and  the  mesh  derivative.  For  this  example,  dqM(£,  q)  is  easily  accessible.  However,  for  2D  and  3D 
problems,  these  maps  are  constructed  using  numerical  algorithms,  and  obtaining  derivatives  of  the 
maps  can  be  very  difficult.  Using  the  preceding  definitions,  we  define  two  methods  for  computing 
the  sensitivity  s(x,q). 

Hybrid-Sensitivity  Equation  Method  (H-SEM) 

Step  1.  Solve  the  transformed  state  equation  (4.1. 15)-(4. 1.16)  for  u>(£,g). 

Step  2.  Solve  the  transformed  sensitivity  equation  (4.1.17)-(4.1.18)  for  s(£,  q). 

Step  3.  Map  s(f,g)  back  to  the  physical  domain  to  obtain  the  sensitivity  by 

s(x,q)  =  s(T(x,q),q). 

Abstract-Semi  Analytical  Method  (A-SAM) 

Step  1.  Solve  the  transformed  state  equation  (4.1.15)-(4.1.16)  for  u>(£,g). 

Step  2.  Solve  (4.1.21)-(4.1.22)  for  the  sensitivity  of  the  transformed  state  p(f,g). 

Step  3.  Map  p(£,g)  back  to  the  physical  domain  to  obtain  the  sensitivity  using 

s(x,q)-p(^,q)-diw(^,q)  d^M{^,q)  dqM(^,q)  . 

Differences  in  the  regularity  of  the  sensitivity  equations  have  already  been  noted.  We  want 
to  examine  the  effect  of  these  differences  on  numerical  approximations  of  the  sensitivities.  The 
following  sections  explore  some  computational  issues  that  are  relevant  to  each  method.  Variational 
formulations  are  described,  and  a  brief  section  outlining  the  discretization  is  presented. 


4.1.6  Variational  Formulations 


We  begin  by  considering  the  transformed  state  equation  in  (4.1.15)-(4.1.16).  Multiplying  by  an 
arbitrary  function  77  G  V  and  integrating  by  parts,  we  have  the  following  integral  equation 

f1  = q2  f  mum  (4.1.26) 

Jo  Jo 

This  equation,  along  with  the  bilinear  form  and  L2-inner  product  defined  in  Section  4.1.2,  produces 
the  variational  form  of  (4.1.15)-(4.1.16).  In  particular,  (4.1.15)-(4.1.16)  is  equivalent  to  the  following 
variational  equation.  Find  u>(-)  €  V  such  that 

a(w,i])  =  q2{f(-,q),T](-))  (4.1.27) 


for  all  q(-)  €  V. 

We  now  turn  to  the  transformed  sensitivity  equation.  Note  that  the  boundary  conditions  (4.1.18) 
are  nonhomogeneous.  In  order  to  simplify  notation,  we  denote  the  right  boundary  condition  of 
(4.1.18)  by  7  =  ^^u)(l).  If  s*(£,  q)  is  defined  by  the  function 

«*(£,?)  =7*.  (4.1.28) 

then  s*(£,q)  €  Hl  (fli )  and  s*  (£,  q)  satisfies  the  boundary  conditions  s*(0)  =  0  and  s*(l)  =  7.  It 
follows  that  v(Z,  q)  —  s(£,  q)  —  s*(£,q)  belongs  to  V  and  solves  the  differential  equation 

~v"(t,q)=  0  (4.1.29) 

with  homogeneous  Dirichlet  boundary  conditions 


v(0,q)  =  0  v(l,g)  =  0.  (4.1.30) 

The  corresponding  variational  equation  is  defined  in  V.  Find  v(-)  €  V  such  that 


a(v,T))=  0 


(4.1.31) 


for  all  tj(-)  €  V.  Once  u(£,  q)  is  computed,  the  transformed  sensitivity  s(£,  q)  is  recovered  using  the 
relationship  s(£,q)  =  v(£,  q)  +  s*(£,q). 

We  now  derive  the  variational  problem  involving  the  sensitivity  of  the  transformed  state,  p(£,q). 
As  noted  in  Section  4.1.5,  the  system  in  (4.1.21)-(4.1.22)  must  be  interpreted  in  the  space  V*.  That 
is,  the  following  equation  holds 


1: 


pwm 


2 q  f  mq)T](m  -  t  S'-iZMOdZ 

Jo  Jo  ’ 

2 qf  ht,q)v(Z)dZ-ri(h 
Jo  q 


(4.1.32) 


for  all  tj(-)  €  V.  Define  the  linear  functional  lq  €  V*  by 


lg(q)  =  2g(/, 77)  -  6i(r})  (4.1.33) 

for  all  rj(-)  6  V.  It  follows  that  equation  (4.1.32)  is  equivalent  to  the  following  variational  problem. 
Find  p(-)  e  V  so  that 


for  all  jj(-)  e  V. 


a(p,ri)  =  lq(r}) 


(4.1.34) 
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Before  turning  to  numerical  issues,  we  note  that  this  example  falls  within  a  general  framework 
that  provides  the  basic  existence,  uniqueness  and  differentiability  of  the  transformed  state  equation. 
Observe  that  the  variational  problem  (4.1.27)  has  the  form 

Aw  =  F(q)  in  V *,  (4.1.35) 

where  A:  V  ->  V*  and  F:  R  V*  are  defined  by 

[Au](v)  =  a(u ,  v)  V  v  e  V  (4.1.36) 

and 

tF(9)](v)=92</(-,  g)M’))  (4.1.37) 

respectively.  If  one  defines  the  function  G :  V  x  R  V*  by 

G(u>,  q)  —  Aw  -  F(g),  (4.1.38) 

then  the  following  result  provides  the  theoretical  foundation  for  the  numerical  methods  presented 
below. 

Theorem  4.1.1.  If  q  £  (1,2),  then  G (w,q)  —  0  has  a  unique  solution  w  =  w(',q).  Moreover, 
p(*,g)  =  9qw(^q)  exists  and  is  the  unique  solution  to  the  sensitivity  equation  (in  V*) 

Ap  -  dqF{q)  —  0.  (4.1.39) 

The  proof  of  Theorem  4.1.1  follows  easily  from  the  implicit  function  theorem.  In  particular,  one 
can  show  that  the  strong  Frechet  derivatives  du,G (w,q)  and  dqG (w,q)  exist  and  are  given  by 

d«,G  (w,q)  =  A,  (4.1.40) 

and 


dqG(w,  q)  =  -dqF(q)  =  -/,(•),  (4.1.41) 

where  /,(•)  is  defined  by  (4.1.33).  Observe  that  the  sensitivity  equation  (4.1.39)  is  equivalent  to  the 
variational  problem  (4.1.34). 

Although  the  framework  defined  by  (4.1.35)  -  (4.1.39)  is  suitable  for  many  elliptic  problems, 
it  is  not  sufficient  for  more  general  shape  sensitivity  problems.  For  example,  the  classical  elliptic 
interface  problem  requires  a  more  general  theory.  We  now  move  to  the  numerical  approximation 
of  the  solutions  to  the  variational  problems.  The  following  section  describes  the  appropriate  finite 
element  spaces  for  approximating  the  transformed  state,  w(£,q),  the  transformed  sensitivity,  s(£,q) 
and  the  sensitivity  of  the  transformed  state,  p(£,q). 

Discretization 

For  the  finite  element  approximation  of  the  variational  equations,  we  begin  by  constructing  the 
grid.  Note  that  the  function  f(£,q)  is  discontinuous  at  the  point  £  =  ~  in  the  computational 
domain.  Hence,  a  grid  point  of  the  mesh  is  placed  at  that  point.  We  partition  the  domain  fix  into 
subintervals  (£/,£j+i)  where  0  =  £0  <  fi  <  <  6  =  \  <  <  ftr+i  =  1.  We  choose 

the  finite  dimensional  subspace  of  V  to  be  the  space  spanned  by  N  piecewise  linear  basis  functions 
denoted 

vN  =  m  e  v  ■.  m  =  f;  (4.1.42) 

3= 1 
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where  hj(£)  is  the  standard  continuous  piecewise  linear  basis  function. 

For  the  computations  presented  in  the  following  sections,  we  use  the  partition  of  flj  developed 
above  for  both  state  and  sensitivity  approximations.  Hence,  the  mesh  1),  for  j  =  0, 

is  used  to  calculate  approximations  to  v(£)  and  p(£,  q).  The  function  u(£)  is  approximated 

in  the  same  manner  as  the  transformed  state.  In  particular,  we  let 

N 

tiN(0  =  ajM0«  (4.1.43) 

j= 1 

N 

PN(Z,Q)  =  '52Tjhj(0  (4.1.44) 

3- 1 


=  (4-1-45) 

i=i 

be  finite  element  approximations  of  u)(£,g),  p(£,g)  and  v(f),  respectively.  Once  vN(£)  has  been 
obtained,  we  compute  the  approximation  of  the  transformed  sensitivity  using 

*2(6  9)  =  ^(e)  +  /(.  (4-1-46) 

where  7 N  is  the  approximation  to  7  =  ^d^w^l^q)  given  by 

~  -~d£WN(l9q).  (4.1.47) 

Note  that  the  subscript  notation  used  on  £$(£,<?)  refers  to  the  use  of  the  H-SEM  approach  for  the 
sensitivity  calculation. 

Remark:  The  use  of  7^  for  7  introduces  error  into  (£,<?)  that  is  independent  of  the  error  in 
^(£>0)*  However,  this  error  significantly  affects  the  accuracy  of  s(f,g),  and  eventually,  that  of 
s(x,q).  It  is  also  noteworthy  that  even  though  the  weak  form  of  (4.1.34)  does  not  require  spatial 
information  about  u)(£),  the  spatial  derivative  d^w(^)  is  required  to  reconstruct  $(£,</)  through 
(4.1.25).  These  issues  play  an  important  role  in  the  numerical  results  presented  in  Section  4.1.7. 

4.1.7  Numerical  Results 

In  this  section,  we  present  numerical  approximations  to  w(x,  q)  and  s(x,q)  using  the  methods  de¬ 
scribed  in  the  previous  section.  We  compare  numerical  approximations  of  s(x,q)  that  result  from 
implementing  each  of  the  methods  discussed  in  Section  4.1.5.  In  particular,  we  compare  results 
obtained  using  H-SEM  with  those  of  the  A-SAM. 

State  Approximations 

All  computations  presented  use  the  same  grids  for  both  sensitivity  approximations.  It  is  important 
to  recall  that  a  node  is  placed  at  £  =  Figure  4.1.1  shows  the  finite  element  approximations 
to  iSAr(f,1.5)  for  various  values  of  N.  These  approximations  converge  rapidly.  Similar  behavior 
is  observed  over  a  range  of  parameter  values.  The  corresponding  approximations  to  wN(x,  1.5), 
obtained  by  transforming  the  finite  element  approximation,  wN(£,q),  back  to  the  physical  domain, 
are  shown  in  Figure  4.1.2.  Comparing  Figures  4.1.1  and  4.1.2,  one  can  see  that  convergence  of  the 
approximations  is  preserved  under  the  transformation  T. 

Figure  4.1.3  shows  the  H^-error  in  w(x,  q)  for  values  of  q  between  1.1  to  1.9.  The  values  of  N 
range  from  3  to  33  and  are  indicated  in  the  legend.  Note  that  the  rate  of  convergence  is  better  for 
q  -»■  1  as  the  quadratic  term  (see  (4.1.19))  in  the  transformed  state  becomes  less  dominant. 
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igure  4.1.1:  Finite  Element  Approximations  to  t/)(£,  1.5) 


Figure  4.1.2:  Approximations  to  w(x ,  1.5) 
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Figure  4.1.3:  i/1-Error  of  wN(x)  for  q  ranging  from  1.1  to  1.9 


Hybrid  Sensitivity  Equation  Method 

In  this  section,  we  present  sensitivity  calculations  obtained  by  applying  the  H-SEM  algorithm. 
Figure  4.1.4  shows  the  convergence  of  the  finite  element  approximations  to  1.5).  Observe  that 
the  entire  error  results  from  the  approximation  of  7  =  ^-d^wil^q)  by  7N  =  :=fd^wN(  l9q).  Since 
the  transformation  T  is  smooth,  the  only  error  in  s^(x,  1.5)  is  due  to  this  approximation.  Figure 
4.1.5  shows  the  convergence  of  s£(x,  1.5)  to  $(#,  1.5). 

Recall  that  wN{£)  is  a  piecewise  linear  approximation.  Thus,  the  finite  element  spatial  derivative 
is  a  piecewise  constant  function.  This  function  is  used  to  approximate  the  spatial  derivative  at  the 
right  boundary  point  f  =  L  Figure  4.1.6  shows  a  piecewise  constant  approximation  used  to  obtain 
an  approximate  boundary  condition  7^.  Hence,  the  error  in  9^u)(l)  results  in  sensitivity  errors  that 
can  be  attributed  to  the  poor  approximation  of  this  boundary  condition.  There  are  techniques  which 
can  be  used  to  obtain  better  approximations  to  the  spatial  derivative  along  the  boundary.  Higher 
order  elements  can  be  used  in  the  transformed  state  calculation,  but  this  can  be  costly  for  2D  and 
3D  problems.  As  an  alternative,  projection  techniques  have  been  developed  to  enhance  the  accuracy 
of  the  spatial  derivative  for  nominal  expense.  We  move  to  the  A-SAM. 

Abstract  Semi- Analytical  Method 

We  turn  our  attention  to  numerical  results  obtained  by  using  the  A-SAM  algorithm  for  sensitivity 
calculations.  First,  note  that  (x)  is  constructed  from  ^(£,<7)  and  dxwN(£)  using  the  relationship. 

*A(z,q)  =PN(^,q)  -  (^)  dzwN(^,q).  (4.1.48) 

Note  that  we  use  the  subscript  A  for  the  sensitivity  approximation  obtained  using  the  A-SAM 
approach.  The  finite  element  approximations  to  jp(£,  1.5)  are  shown  in  Figure  4.1.7  for  various 
values  of  N .  Since  the  sensitivity  equations  are  linear,  the  approximations  pN{£,  1.5)  converge  as 
expected.  When  constructing  1.5)  from  (4.1.48),  the  piecewise  constant  approximation  of 
dxwN(£)  produces  discontinuities  in  s^(x)  as  shown  in  Figure  (4.1.8).  These  discontinuities  occur 
at  points  of  the  physical  domain  which  correspond  to  mesh  nodes  of  the  computational  domain  lying 
in  the  interval  [^,  1).  Note  that  the  expressions  for  the  mesh  derivatives  in  (4.1.48)  are  “hard- wired” , 
continuous  functions,  and  the  finite  element  approximations  to  p(£,  q)  are  continuous.  It  follows  that 
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Figure  4.1.7:  Finite  Element  Approximations  to  p(£,  1.5) 


||s(x,g)  -  s^(x,q)\\L2  0  as  N  ->•  oo.  However,  one  does  not  get  convergence  in  the  energy  (Hq) 

norm  since  dxs^(x,q)  is  not  in  V. 

4.1.8  Conclusion 

Each  of  the  methods  described  above  present  computational  challenges.  The  Hybrid  SEM  requires 
gradient  information  along  the  boundary  of  the  computational  domain.  Even  in  the  ID  example,  this 
information  becomes  critical  to  accurate  sensitivity  calculations.  Approximating  gradients  along  the 
boundary  only  gets  more  challenging  in  2D  and  3D  problems.  However,  the  Abstract  SAM  requires 
accurate  gradient  information  within  the  computational  domain.  This  example  clearly  illustrates 
the  serious  contamination  of  sensitivity  approximations  that  can  occur  if  the  approximate  gradients 
are  inaccurate  or  axe  not  sufficiently  smooth.  Moreover,  the  issue  of  calculating  derivatives  of  mesh 
maps  is  not  addressed  in  this  ID  example.  Those  were  analytically  computed  and  “hard-wired”  into 
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Figure  4.1.9:  Hl  Errors  for  Sensitivity  Calculations 
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the  computations.  In  practical  problems,  calculating  these  derivatives  can  be  an  extremely  difficult 
task  as  well  as  a  source  of  computational  error.  Both  methods  require  accurate  gradient  information 
from  the  transformed  state.  In  the  case  of  the  Hybrid  method,  this  is  the  only  major  stumbling 
block  to  obtaining  reliable  sensitivity  calculations.  In  contrast,  accurate  gradient  information  may 
not  be  sufficient  to  obtain  good  sensitivity  approximations  using  the  Abstract  SAM  method  as  the 
need  for  accurate  mesh  derivatives  may  overshadow  the  entire  process. 
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4.2  Fast  Algorithms  for  Computing  Functional  Gains 

Functional  gains  are  kernels  of  feedback  operators  that  result  from  distributed  parameter  control 
problems  defined  by  partial  and  functional  differential  equations.  These  gains  offer  insight  into 
issues  such  as  sensor /actuator  placement  and  controller  reduction.  To  be  practical,  one  must  be 
able  to  compute  these  kernels  for  a  wide  variety  of  partial  differential  equations  in  2  and  3  spatial 
dimensions.  Standard  “approximate-then-design”  approaches  have  been  very  useful  for  small  ID 
problems.  However,  in  2D  and  3D  problems,  the  size  of  the  approximating  systems  limits  this 
method  as  a  practical  computational  tool.  Therefore,  alternative  methods  are  needed. 

As  a  first  step  in  the  direction  we  have  investigated  the  possibility  of  using  “direct”  approxi¬ 
mations  of  Riccati  and  Chandrasekhar  partial  differential  equations  that  define  the  kernels.  In  this 
short  note  we  illustrate  this  idea  by  using  finite  element  schemes  to  solve  these  partial  differential 
equations.  One  important  observation  about  this  direct  approach  is  that  it  allows  for  the  possi¬ 
bility  of  using  parallel  and  adaptive  computational  tools.  In  addition,  we  show  ^hat  modifications 
of  “standard”  finite  element  schemes  can  often  improve  the  speed  and  accuracy  of  the  old  indirect 
schemes.  The  goal  of  the  current  note  is  to  illustrate  the  basic  ideas  and  provide  some  introduction 
to  the  theoretical  and  computational  issues  that  arise  in  this  problem  area. 


4.2.1  Dirichlet  Boundary  Control 

The  one-dimensional  heat  equation  with  Dirichlet  boundary  controls  is  defined  by  the  initial  bound¬ 
ary  value  problem: 


zt(t,x)  =  ezxx(t,x), 
z(t,  0)  =  u0(t),  z(t,  1)  =  ui(t),  > 
z(0,x)  =  z0(x), 


(4.2.49) 


where  e  is  the  thermal  diffusivity  of  the  material. 

In  order  to  take  advantage  of  distributed  parameter  control  theory,  we  first  formulate  the  bound¬ 
ary  value  problem  (4.2.49)  as  an  abstract  state  space  model.  In  particular,  we  write  the  problem  as 
a  system  of  the  form 


z(t)  =  Az{t)  +  Bu(t), 
z(0)  =  z0, 

where  A  and  B  are  defined  below.  Let  A  be  the  differential  operator 

[M]{X)  =  C0(a:)’ 


(4.2.50) 


(4.2.51) 


defined  on  the  domain  V{A)  =  Hq(0,  1)  n  H2{ 0, 1),  and  let  D  denote  the  Dirichlet  map  D  :  K2  -4 
L2( 0, 1)  given  by 


[Dw](x)  =  (1  —  x)uo  +  xui.  (4.2.52) 

The  operator  A  :  L2( 0, 1)  — >■  [D(A*)}'  is  the  lifting  of  A  to  L2( 0, 1)  defined  by 

I4/KV0  =  =  e  [  f(x)ip"(x)dx,  (4.2.53) 

Jo 

for  each  ip  €  T>(A*)  =  Hq  (0, 1)  fl H2(0, 1).  Finally,  the  input  operator  B  :  R2  -»•  [D(A*)]'  is  given  by 

B  =  -AD.  (4.2.54) 
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The  adjoint  of  B  is  usually  easier  to  represent  than  the  operator  itself;  moreover,  it  is  the  adjoint 
that  one  needs  in  finite  element  approximations.  In  particular,  one  has 


(Bu,  (f>) 


—e  /  [(1  -  x)u0  +  xui](j)n (x)dx 
Jo 

1  rl 

-e[(l  -  x)u0  +  xui]<ff(x)\0  +  e(ui  -  u0)  /  </>'(x)dx 

Jo 

-eu^iD  +  eGuo^iO) 

(u,B*<j>)u  2. 


Hence,  B*  :  V(A*)  -4  M2  is  given  by  B*<j>  —  0),  1)]T. 

We  consider  a  weighted  LQR  control  problem  defined  by  minimizing  the  cost  function 


m-)) 


-f 


{(y(t),y(t))  +  {u{t),u(t))}e  “  dt 


(4.2.55) 


subject  to  the  constraint  (4.2.50),  where  y(t)  =  Cz(t)  is  a  controlled  output  function.  It  is  well  known 
(under  suitable  conditions  on  C)  that  the  optimal  control  for  this  problem  is  given  by  feedback  of 
the  form 


u{t)  =  Kz(t)  =  - B*Uz(t ),  (4.2.56) 

where  K  is  a  feedback  gain  operator.  Here,  II  is  the  (weak)  solution  to  the  algebraic  Riccati  equation 
(ARE) 


(A  +  aiy n  +  n(i  +  ai)  -  hbb*u  +  q  =  o,  (4.2.57) 


and  Q  =  C*C. 

It  is  possible  to  show  that  the  feedback  operator  K  maps  L2( 0, 1)  to  R2  and  has  the  form 

Kcj>  -  (4.2.58) 

where  the  kernel  h(-)  belongs  to  L2(0, 1).  This  kernel  is  called  the  functional  gain  and  is  useful 
in  many  applications.  In  previous  research,  Burns  and  Rubio  use  functional  gains  to  optimally 
place  sensors.  Similar  ideas  have  been  used  to  construct  low-order  observers.  A  first  step  in  both 
applications  is  the  development  of  fast  and  accurate  numerical  schemes  for  computing  h(-).  This  is 
the  central  theme  of  this  work. 

4.2,2  Riccati  Methods 

Conceptually,  there  are  several  ways  to  compute  the  functional  gain.  One  basic  approach  is  to 
approximate  the  ARE  (4.2.57)  and  obtain  an  approximate  Riccati  operator  n^.  The  approximate 
feedback  operator  is  given  by 


Kn  =  -[B*]NnN, 


(4.2.59) 


and  this  feedback  law  yields  an  approximation  hN(-)  to  the  functional  gain  h(-). 

A  more  direct  approach  is  to  first  obtain  a  representation  for  the  Riccati  operator  of  the  form 


[n <f>](x)  =  f  p(x,0<t>( 0d£, 
Jo 


(4.2.60) 
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and  to  then  solve  for  p(a;,f).  For  certain  systems  (distributed  control  and  Neumann  boundary 
control),  Lions  uses  the  Schwartz  Kernel  Theorem  to  obtain  a  representation  of  the  form  (4.2.60). 
Lions’  result  does  not  apply  directly  to  the  Dirichlet  boundary  control  problem.  However,  King 
extended  this  result  to  the  one-dimensional  case  considered  here. 

Theorem  4.2.1.  Assume  Q  :  L2{ 0, 1)  L2{ 0, 1)  is  bounded .  If  II  =  II*  is  the  unique ,  nonnegative- 
definite  solution  to  the  ARE  (4.2.57),  then  n  is  Hilbert-Schmidt.  Moreover,  there  exists  a  function 
p(x,£)  such  that  n  has  the  representation  (4.2.60)  where  the  kernel  p{x,()  satisfies  the  condition : 

p(x,0=p(t,x)  e  C^MxtO.l]).  (4.2.61) 

In  order, to  apply  this  result,  we  assume  that  C  :  L2( 0, 1)  R  has  the  form 

C<t>( *)  =  f  c(x)(f>(x)dx ,  (4.2.62) 

Jo 

where  c(x)  E  L2( 0, 1).  Note  that  the  operator  Q  =  C*C  has  the  form 

[Q(f>](x)  =  c(x)  f  c(f)<£(f)d£ 

Jo 

=  / 

Jo 

where  q(x,£)  is  the  kernel  for  Q. 

The  following  theorem  extends  the  results  in  Lions  to  the  Dirichlet  boundary  control  problem 
considered  here. 


Theorem  4.2.2.  The  kernel  p(x,£)  of  the  operator  n  is  a  weak  solution  to  the  Riccati  partial 
differential  equation  (R-PDE) 


(Ax  +  A/:  +  2al)p(x,£)  -  -^p(0,Z)^p(x,0) 

+  +  =  o. 

In  particular,  p(-,  •)  €  Hq  ([0, 1]  x  [0, 1])  satisfies  the  boundary  conditions 

P(0,0  =  J»(1,0  =  0,  and 
p(x,0)  =  p(x,  1)  =  0, 

and  the  symmetry  condition 

p(x,4)  =  p((,x),  Vx,{e(o,i). 


(4.2.63) 


(4.2.64) 


(4.2.65) 


Proof.  The  proof  follows  the  same  basic  idea  used  by  Lions  for  the  case  of  distributed  control 
with  zero  Dirichlet  boundary  conditions.  Property  (4.2.65)  follows  from  II*  =  II;  that  is, 

p*]oo  =  fpMmds  =  f  p(t,x)mds  =  pr*](*). 

Jo  Jo 

From  Theorem  3.1,  p(-,  •)  is  in  C'1([0, 1]  x  [0, 1]).  Suppose  that,  without  loss  of  generality,  p( 0,  fo)  >  0 
for  some  fixed  value  of  £o  €  (0, 1).  Since  p(0,  •)  €  C'1(0, 1),  there  exists  a  neighborhood  No  C  (0, 1) 
containing  f0  such  that  p(0,£)  >  0  for  all  £  €  N0.  It  follows  that  II 0  €  V(A)  and  hence  [110] (0)  =  0 
for  all  0  E  i2(0, 1).  If  <f>(x)  —  1  on  No  and  <f>(x)  =  0  elsewhere,  then 

[n0](o)  =  fP( 0,00m  =  /  p(o,0^  #  0, 

Jo  Jn0 
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which  contradicts  the  equality  [II<£](0)  =  0.  Therefore,  p(0,  £)  =  0  for  all  £  €  (0,1).  The  other 
boundary  conditions  in  (4.2.64)  are  established  by  using  analogous  arguments. 

To  obtain  (4.2.63),  we  substitute  the  representation  of  the  Riccati  operator  into  (4.2.57).  Recall 
that  the  weak  form  of  the  ARE  (4.2.57)  is  defined  by 

(U<f>,A*rp)  +  (A*<i>,nil>)+2a(n<i>,  ifi)  -  (B*U(f>,  B*Hil))  +  {C<t>,Cil>)  =  0,  (4.2.66) 

for  each  <j>,ip  €  V(A).  In  particular,  the  first  term  in  (4.2.66)  is  given  by 

(n<M»  =  j f1  nx)dx 

Let  V  =  JTo([0, 1]  x  [0, 1])  and  define  Ax  :  V  -4  V*,  by 

[A.p(*,  {)](»>(•,•))  =  ~  JQ  f0  §^P(x’0-^v(x,Odxd£. 

The  operator  A$  is  defined  similarly.  The  term  (Cp,  Cip)  takes  the  form 

(ere*,*)  =  m,ip) 

~  J0\J0  Q(x’0<Km  ip(x)dx 

=  f  [  q(x,04>((;)ip(x)dxd£. 

Jo  Jo 

Finally,  since  p(-,  •)  €  C1([0, 1]  x  [0, 1]),  it  follows  that 
(UBB^tp)  =  (B*n<t>,B*n*il>)R2 

=  [B*  £  p(x,  &#€)#]  T  [**  fo  P(*’  0^(x)dx] 


=  So  So  [kP(0»O^p(«.0) 

+^P(l»0^(*.l)]  4>(Oi>(x)dxd£. 

Therefore,  p(-,  •)  €  V  =  Hq  ([0, 1]  x  [0, 1])  satisfies  the  variational  problem:  Find  p(-,  •)  G  V  such  that 
for  all  v(-,  -)  €  V 

n1  d  d  f1  f1  d  d 

fap(X’®fa;V(x’®dxdZ-J0  J  Q£P(x’Z)Q£V(x’Odxd£ 

+  2a  [  [  p(x,t)v(x,£)dxd£ 

i  i  !  (4.2.67) 

~i  L  [|(0’ol(x,0)+f(1’e)i(:r’1)]u(a:’e)da:^ 

+  f  q(x,Z Mx,0  =  o. 

Jo 

Thus,  the  kernel  p(x,£)  is  a  weak  solution  to  (4.2.57).  □ 

The  feedback  gain  operator  has  the  form 


K(j>  = 


mm 
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Moreover,  it  has  been  shown  that  B*  can  be  moved  inside  the  integral  so  that 

K4>  =  f1 -[B*xP](x,0<Km- 
Jo 

Thus,  the  functional  gain  h(£)  is  given  by 

MO  =  -[£>](*,  0  =  [  dp( 0, Q/dx,  -dp(  1,  Q/dx  ]T. 


4.2.3  Chandrasekhar  Methods 


Chandrasekhar  methods  have  been  used  to  bypass  the  step  of  solving  for  the  Riccati  operator. 
Formally,  the  Chandrasekhar  equations  are  defined  by 


K{t)  =  B*L*(t)L{t), 

L(t )  =  -L(t)[A  +  al  -  BK(t)], 


(4.2.68) 


with  initial  conditions 


K(T)  =  0  and  L(T)  =  C.  (4.2.69) 

Chandrasekhar  equations  have  been  applied  to  a  number  of  distributed  parameter  systems.  Under 
suitable  conditions,  one  can  show  that  K(t)  K  ast  -oo.  For  the  problem  considered  here,  the 
Riesz  Representation  Theorem  yields  a  functional  gain  h(-)  so  that 


[K(p\  =  f  h{x)4>{x)dx.  (4.2.70) 

Jo 

One  can  discretize  the  Chandrasekhar  equations  in  the  same  manner  as  the  discretization  of 
the  Riccati  equations.  In  particular,  one  can  construct  approximations  of  the  system  and  solve  the 
corresponding  ordinary  differential  equations 


KN(t)  =  [£*]"[£*]"(*)£"(<), 

LN(t)  =  -LN(t)[AN  +oJ-  BNKN(t)}, 


(CEN) 


with  initial  conditions 


Kn(T)  =  0  and  LN(T)  =  CN . 


(4.2.71) 


This  approach  has  been  used  by  a  number  of  researchers.  However,  one  loses  spatial  information 
that  can  be  helpful  in  developing  adaptive  finite  element  methods  for  direct  computation  of  h{-).  A 
second  approach  is  to  derive  the  Chandrasekhar  partial  differential  equations  similar  in  spirit  to  the 
Riccati  partial  differential  equation  (4.2.63). 

One  formulates  this  system  of  partial  differential  equations  for  the  functional  gain  by  first  as¬ 
suming  representations  of  the  operators  in  the  Chandrasekhar  equations  (4.2.68)  and  then  proving 
the  existence  of  these  kernels.  In  particular,  we  seek  kernels  h(t,£)  and  l(t,  £)  such  that 

K(t)4>  =  [  h(t,  mm  (4.2.72) 

Jo 

and 

L(t)ct>  =  /  (4.2.73) 

Jo 
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respectively. 

The  kernels  ft(-,-)  and  /(*,*)  satisfy  the  following  Chandrasekhar  partial  differential  equations 
(C-PDE) 


dh 

dt 


m 

at 


m  =  -w+anmo+hfrwijfro,  J 


with  initial  conditions 


(4.2.74) 


h(T,x)  =  0  and  l(T,x)  =  c(x),  (4.2.75) 

where  c(x)  is  the  kernel  of  the  operator  C .  Also,  h(-,  •)  and  /(•,  •)  satisfy  the  boundary  conditions 

h(t,0)  =  h(t,  1)  =  0  (4.2.76) 


and 


l(t,0)  =  l(t,l)  =  0,  (4.2.77) 

respectively. 

Note  that  the  term  B*l  produced  a  boundary  value  so  that  (4.2.74)-(4.2.77)  is  a  rather  complex 
boundary  value  problem.  In  particular,  the  boundary  terms  B*  acting  on  the  kernel  l(t,  •)  yields 

pma  =  [j^M),  -^i(Mi]  • 

The  direct  method  involves  solving  (4.2.74)  for  /i(£,  ■)  and  l(t ,  •),  and  using  the  fact  that  as  t  — oo, 
h(t,£)  -+  MO  (the  desired  functional  gain).  To  complete  this  program,  one  needs  to  establish 
existence  and  regularity  of  h(t,  f)  and  /(£,£).  In  order  to  limit  the  length  of  this  paper  we  do  not 
address  these  issues  here. 


4.2.4  Appr  oximat  ions 

We  compute  the  functional  gains  by  applying  several  finite  element  schemes  to  both  approaches 
discussed  above.  In  one  case,  we  approximate  the  system  (AN  ,BN ,QN ,RN)  and  solve  the  finite 
dimensional  Riccati  and  Chandrasekhar  equations.  We  also  present  results  for  direct  finite  element 
approximation  of  the  Chandrasekhar  partial  differential  equations. 

Once  a  finite  element  subspace  is  selected,  the  first  approach  is  to  build  matrix  representa¬ 
tions  of  the  operators  of  the  system  (4.2.50)  with  cost  function  (4.2.55).  In  particular,  we  obtain 
^  «  A,  Bn  «  J3,  CN  «  C  and  RN  =  R  =  L  These  matrix  representations  yield  approximate 
Riccati  equations  (or  Chandrasekhar  equations)  which  can  be  solved  using  an  appropriate  numerical 
algorithm.  The  approximate  solution  11^  or  KN  yields  an  approximation  hN( •)  of  h(-). 

We  consider  three  finite  element  methods,  and  all  methods  employ  piecewise  linear  elements.  It 
is  important  to  review  these  methods  because  the  three  schemes  differ  only  in  the  implementation 
of  the  boundary  conditions. 

Recall,  the  global  basis  functions  are  defined  by  the  piecewise  linear  “hat”  functions 

{(x  —  Xi-i)/h,  1  <  X  <  Xi 

(Xi+1  -  x)/h,  Xi  <  X  <  Xi+ 1  (4.2.78) 

0,  otherwise , 

,  N ,  and  AT  +  1  is  the  number  of  subintervals  of  (0, 1).  At  x  =  0 
given  by 

f  (h  —  x)//i,  0  <  x  <  h 

\  0,  otherwise , 


hi(x)  = 

where  h  =  i/(N  4*  1)  for  i  =  1, . . . 
and  x  —  1,  ho(x)  and  h^+i(x)  are 

h0(x) 
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and 


W(»)  =  { 


XN  <  X  <  l 

otherwise , 


respectively. 

The  first  method  is  “standard  finite  elements.”  The  second  method  was  used  by  Burns  and  Kang 
in  their  study  of  Burgers’  equation  and  uses  all  N  4-  2  basis  functions.  The  third  method  is  based 
on  Nitsche’s  approximation  as  discussed  in  Section  10.1  of  Lasiecka  and  Triggiani. 

The  standard  finite  element  method  is  the  most  commonly  discussed  scheme.  We  find  the  solution 
to  the  problem  in  the  subspace  Vq  C  Hq  (0, 1),  where  Vq  is  the  space  spanned  by  hi{x),  i  —  1, . . .  ,N. 

In  Vq,  the  matrix  representation  for  A  is  Aff  =  (M^)~1Sq  where  Mq  and  Sq  are  the  mass 
matrix  and  stiffness  matrix,  respectively.  These  matrices  are  given  by 


'4  1  0  •••  0‘ 

'  -2  1  0  •• 

•  0  ' 

1  4  1  •••  0 

1  -2  1  •• 

•  0 

0  1  4  •••  0 

oN  _  6 

’  ^  "  h 

0  1  -2  •• 

•  0 

0  0  0  •••  4 

0  0  0  •• 

•  “2  . 

The  approximation  Bq  has  the  representation  ( MN )  1Bq  where 


“10“ 
0  0 

0  0 
0  1 


(4.2.79) 


The  remaining  finite  element  schemes  attempt  to  adjust  for  the  fact  that  the  solutions  are  not 
zero  on  the  boundaries.  Bramble,  et.  al.,  consider  these  kinds  of  subspaces  because  they  give  optimal 
convergence  rates  for  nonhomogeneous  boundary  conditions.  Burns  and  Kang  employ  all  iV+2  basis 
functions  and  hence  project  the  problem  onto  the  space  Vh  spanned  by  hi(x ),  i  =  0, 1, . . .  ,N  +  1. 

The  resulting  matrix  representations  for  AN  and  BN  are  of  size  (N  +  2)  x  (N  +  2)  and  (N  + 
2)  x  2,  respectively.  Thus,  AN  has  the  representation  ( MN)~lSN  and  BN  has  the  representation 
{MN)~lBN  where 


'210 

...  o  ' 

O 

O 

o 

H 

1 

1  4  1 

...  o 

1  -2  1  •••  0 

0  14 

...  0 

qN  _  e 

’  5  h 

0  1  -2  •••  0 

,000 

...  2 

0  0  0  •••  -1 

and  BN  is  the  (N  +  2)  x  2  analog  of  (4.2.79). 

Lasiecka  and  TViggiani  suggested  that  Nitsche’s  approximation  would  be  useful  in  the  LQR 
problem  because  of  improved  accuracy  of  the  Dirichlet  boundary  conditions.  This  scheme  produces 
optimal  convergence  rates  in  the  solution  of  the  open-loop  system.  As  we  shall  see  below,  the  Nitsche 
scheme  produces  excellent  approximations  to  the  functional  gains.  The  Nitsche  method  begins  with 
the  finite  element  scheme  used  by  Burns  and  Kang.  Nitsche’s  approximation  A@  of  the  A  operator 
is  defined  by 

(Ap  zN ,  <j>N)  =  a(zN,<t>N)-(^zN,<t,N)r-(zN,^<pN)r-fih-1(zN^N)r,  (4.2.80) 
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where  (3  >  0  is  sufficiently  large,  T  is  the  boundary  of  the  domain  and  a(-,  •)  is  the  bilinear  form 
such  that 


a(0,V>)  =  —  f  <j>f(x)'i/j,(x)dx. 

Jo 


Rearranging  the  terms  of  Nitsche’s  approximation  A$  and  using  the  fact  that  2^(0)  =  uo  and 
zN(  1)  =  tii,  the  only  matrices  that  are  changed  are  the  stiffness  and  the  control  matrices.  The 
resulting  stiffness  matrix  is  the  (N  +  2)  x  (N  +  2)  matrix 


(4.2.81) 


-1-13 

0 

0  ••• 

0  ' 

1 

-2 

1  ••• 

0 

qN  _  € 

SP  ~  h 

0 

1 

-2  ... 

0 

l 

o  • 

0 

•  o 

-1-/3  . 

and  the  control  matrix  is  given  by 


(4.2.82) 


Again,  Q  is  approximated  by  projecting  onto  Vh  producing  the  approximate  ARE  defined  by 

{An  +  aIN)*UN  +  Un(An  +  aIN)-nNBN{BN)*IlN  +  Qn  =  0,  (4.2.83) 

It  is  clear  that  one  can  generate  finite  element  approximations  to  the  Riccati  and  Chandrasekhar 
partial  differential  equations  by  projecting  these  systems  onto  the  various  finite  element  spaces.  In 
order  to  conserve  space,  we  focus  on  the  Chandrasekhar  PDE  (4.2.74)-(4.2.77).  We  turn  now  to 
numerical  results  that  compare  the  schemes  and  methods  mentioned  in  this  section. 


4.2.5  Numerical  Results 

In  this  section  we  present  numerical  comparisons  for  the  methods  described  above.  Two  examples 
are  presented.  The  first  example  compares  the  results  obtained  by  applying  the  indirect  method.  We 
approximate  the  functional  gains  using  the  three  finite  element  schemes  mentioned  in  the  previous 
section  via  the  approximate  ARE  (4.2.83). 

The  second  example  focuses  on  the  Chandrasekhar  PDE  system.  We  compare  the  standard  finite 
element  scheme  with  Nitsche’s  method.  In  particular,  we  use  the  standard  finite  element  scheme  to 
solve  the  Chandrasekhar  PDE’s  (4.2.74)-(4.2.77),  and  compare  this  result  to  the  Nitsche  solution  of 
the  ARE  (4.2.57). 

Example  1.  In  this  example,  we  take  Q  —  I  in  (4.2.55),  and  set  e  =  1/60,  a  =  0.3.  Figures 
4.2.10  find  4.2.11  show  plots  of  the  approximate  functional  gains  for  all  three  methods.  Figure  4.2.10 
corresponds  to  control  at  x  =  0  and  Figure  4.2.11  contains  plots  for  the  functional  gains  for  control 
at  x  =  1.  Observe  that  at  N  =  8,  the  Nitsche  scheme  “corrects”  the  errors  at  the  boundary  that 
result  from  a  direct  application  of  the  Burns-Kang  scheme.  The  “converged”  solution  is  represented 
by  the  N  =  128  standard  finite  element  solution.  The  closed-loop  response  is  plotted  in  Figure 
4.2.12.  Figure  4.2.13  provides  comparison  of  the  closed-loop  trajectories  at  time  t  =  10  seconds. 
Note  that  the  Nitsche’s  scheme  provides  more  accurate  responses  than  the  other  methods.  This 
observation  has  potential  applications  to  controller  reduction  schemes. 

Example  2.  In  this  example,  c(x)  is  defined  by 
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0.6  <  £  <  0.8 

otherwise , 


and  e  =  0.1,  a  —  0.  The  standard  finite  element  scheme  is  applied  to  the  Chandrasekhar  PDE 
system  (4.2.74)-(4.2.77).  Nitsche’s  scheme  is  used  to  compute  hN(x)  via  the  ARE  (4.2.83).  Figures 
4.2.14  and  4.2.15  show  plots  of  the  approximate  functional  gains  for  both  schemes.  The  “con¬ 
verged”  solution  is  represented  by  the  same  scheme  used  in  the  previous  example.  Note  that  both 
schemes  produce  excellent  approximations  to  the  functional  gains.  It  is  important  to  note  that  the 
convergence  in  space  ( N  oo)  and  time  ( t  -oo)  of  hN(t,£)  to  /i(£)  is  “fast  enough”  so  that 
approximations  of  the  Chandrasekhar  PDE’s  are  competitive  with  the  other  indirect  approaches. 


Figure  4.2.10:  Functional  Gain  at  x  =  0,  N  =  8. 


Figure  4.2.11:  Functional  Gain  at  x  —  1,  N  —  8. 
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Figure  4.2.12:  Solution  of  Closed-Loop  System,  t  =  0  to  t  =  20. 

4.2.6  Conclusion  and  Future  Research 

Although  the  numerical  results  presented  here  are  not  sufficient  to  provide  conclusive  information 
on  the  most  efficient  method  for  approximating  functional  gains,  they  do  provide  insight  into  these 
matters  and  they  demonstrate  the  feasibility  of  the  various  schemes. 

It  is  interesting  to  note  that  Nitsche’s  approximation  produces  excellent  approximations  to  the 
functional  gains.  Hence,  accuracy  is  improved  by  use  of  the  Nitsche  scheme.  However,  the  Nitsche 
scheme  was  considerably  slower  than  the  standard  finite  element  method.  This  was  especially  notice¬ 
able  as  N  got  larger.  This  is  an  issue  that  should  be  explored.  In  addition,  the  improved  accuracy 
obtained  by  Nitsche’s  scheme  was  not  uniform  in  N.  In  some  cases,  the  standard  finite  element 
scheme  was  better  at  small  N.  Consequently,  the  tradeoffs  between  speed  and  accuracy  are  not  yet 
settled. 

Indirect  methods  (approximate  first)  have  been  useful  in  dealing  with  ID  and  2D  problems. 
However,  these  indirect  methods  do  not  take  advantage  of  the  structure  of  the  problem  and  do 
not  allow  for  adaptive  finite  element  refinements.  The  Riccati  and  Chandrasekhar  PDE’s  offer  the 
potential  for  improvements  in  these  areas.  These  issues  need  to  be  explored. 
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(x)nM 


Figure  4.2.14:  Functional  Gain  at  x  =  0,  N  =  8. 
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Figure  4.2.15:  Functional  Gain  at  x  =  1,  N  =  8 


4.3  Significance  and  Potential  Application  to  Air  Force  Prob¬ 
lems 

Hybrid  Transformation  Methods  for  Sensitivity  Computations  -  Accurate  sensitivity  calcu¬ 
lations  play  an  important  role  in  the  analysis  and  optimization  of  engineering  systems.  Sensitivities 
can  be  used  to  compute  gradients  in  optimization-based  design.  We  focused  on  variational  methods 
for  computing  state  sensitivities.  These  schemes  make  use  of  sensitivity  equation  methods.  However, 
there  are  a  variety  of  ways  to  implement  sensitivity  equation  methods,  and  these  variations  yield 
algorithms  with  different  convergence  properties.  We  considered  two  specific  methods.  The  first  is 
based  on  transforming  the  state  equation  to  a  fixed  computational  domain  and  then  deriving  its 
sensitivity  equation.  Once  the  state  and  sensitivity  systems  are  solved,  the  solutions  are  mapped 
back  to  the  physical  domain.  The  second  approach  transforms  both  state  and  sensitivity  equations, 
solves  the  transformed  equations  and  maps  these  solutions  back  to  the  physical  domain.  There  are 
benefits  and  drawbacks  to  each  method.  Indeed,  it  is  not  always  obvious  which  scheme  is  best  for  a 
given  problem.  Many  questions  needed  to  be  addressed  before  a  complete  theory  could  be  developed. 

Each  of  the  methods  present  computational  challenges.  The  Hybrid  SEM  requires  gradient  in¬ 
formation  along  the  boundary  of  the  computational  domain.  This  information  becomes  critical  to 
accurate  sensitivity  calculations.  However,  other  SEM  often  require  accurate  gradient  information 
within  the  computational  domain.  This  work  clearly  illustrated  that  serious  contamination  of  sen¬ 
sitivity  approximations  can  occur  if  the  approximate  spatial  gradients  are  inaccurate  or  are  not 
sufficiently  smooth.  Moreover,  the  Hybrid  SEM  avoids  the  need  to  calculate  derivatives  of  mesh 
maps.  In  the  case  of  the  Hybrid  method,  accurate  gradients  is  the  only  major  stumbling  block  to 
obtaining  reliable  sensitivity  calculations.  In  contrast,  accurate  gradient  information  may  not  be 
sufficient  to  obtain  good  sensitivity  approximations  using  the  other  SAM  methods,  since  the  need 
for  accurate  mesh  derivatives  may  overshadow  the  entire  process.  This  work  provided,  for  the  first 
time,  an  understanding  of  the  role  that  spatial  smoothness  plays  in  developing  accurate  numerical 
methods  for  sensitivity  computations. 

These  results  will  play  an  important  role  in  developing  fast  and  accurate  sensitivity  analysis  tools 
with  Air  Force  applications  such  as  wing-body  design  optimization,  combustion  control  and  COIL 
lasers. 

Fast  Algorithms  for  Computing  Functional  Gains  -  Functional  gains  are  kernels  of  feedback 
operators  that  result  from  distributed  parameter  control  problems  defined  by  partial  and  functional 
differential  equations.  These  gains  offer  insight  into  issues  such  as  sensor/actuator  placement  and 
controller  reduction.  To  be  practical,  one  must  be  able  to  compute  these  kernels  for  a  wide  variety 
of  partial  differential  equations  in  2  and  3  spatial  dimensions.  Standard  “approximate-then-design” 
approaches  have  been  very  useful  for  small  ID  problems.  However,  in  2D  and  3D  problems,  the 
size  of  the  approximating  systems  limits  this  method  as  a  practical  computational  tool.  Therefore, 
alternative  methods  were  needed. 

As  a  first  step  in  the  direction  we  developed  “direct”  approximations  of  Riccati  and  Chan¬ 
drasekhar  partial  differential  equations  that  define  these  kernels.  This  direct  approach  allows  for 
the  possibility  of  using  parallel  and  adaptive  computational  tools.  In  addition,  we  established  that 
by  modifying  “standard”  finite  element  schemes  one  can  improve  the  speed  and  accuracy  of  the  old 
indirect  schemes.  In  particular,  so  called  “Nitsche  methods”  produces  greatly  improved  approxi¬ 
mations  to  the  functional  gains.  Although  the  standard  Nitsche  scheme  improved  accuracy,  it  was 
considerably  slower  than  other  finite  element  methods.  This  was  especially  noticeable  as  the  number 
of  elements  grew.  Riccati  and  Chandrasekhar  PDE’s  offers  the  potential  for  improvements  in  this 
areas. 

These  methods  have  already  proven  to  be  useful  in  the  control  of  2D  thermal  processes  and  when 
extended  to  Navier-Stokes  equations  will  provide  the  basic  tools  needed  to  accomplish  practical 
sensor/actuator  placement  and  controller  reduction  for  flow  control.  In  addition,  these  tools  will  be 
useful  in  several  other  DOD  applications;  including  the  design  and  control  of  materials  manufacturing 
processes  and  control  of  combustion. 
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Chapter  7 

Interactions  and  Transitions 


One  of  the  major  components  of  this  effort  is  active  cooperation  and  coordination  with  Air  Force 
laboratories  and  with  industrial  partners.  The  ASSERT  students  are  active  participants  in  the 
PRET  Center  and  frequently  interact  with  our  industrial  partners.  In  addition,  these  students  have 
been  invited  to  present  their  research  at  national  and  international  conferences. 


7.1  Participation  and  Presentations  at  Meetings 

Kevin  Hulsing 

1.  Sipe  Consolidation  in  Tire  Design ,  Industrial  Mathematics  Modeling  Workshop  for  Graduate 
Students,  N.C.  State  University,  Raleigh,  NC,  August  1997. 

2.  Methods  for  Computing  Feedback  Operators  for  LQR  Control ,  University  of  Trier,  Trier,  Ger¬ 
many,  February  1998. 

3.  Methods  for  Computing  Feedback  Operators  for  LQR  Control ,  Fourth  SIAM  Conference  on 
Control  &  Its  Application,  Jacksonville,  FL,  May  1998. 

4.  Methods  for  Computing  Feedback  Operators  for  LQR  Control ,  Computation  and  Control  VI, 
Montana  State  University,  Bozeman,  MT,  August  1998. 

Lisa  Stanley 

1.  Two  Projection  Methods  for  Approximating  Derivatives  of  Computed  Solutions  to  PDEs,  Uni¬ 
versity  of  Trier,  Trier,  Germany,  March  1998. 

2.  Two  Projection  Methods  for  Approximating  Derivatives  of  Computed  Solutions  to  PDEs,  Fourth 
SIAM  Conference  on  Control  &  Its  Applications,  Jacksonville,  FL,  May  1998. 

3.  A  Comparison  of  Two  Methods  for  Computing  Sensitivities ,  Control  and  Computation  VI, 
Montana  State  University,  Bozeman,  MT,  August,  1998. 

4.  A  Comparison  of  Two  Methods  for  Computing  Sensitivities ,  Sixth  SIAM  Conference  on  Opti¬ 
mization,  Atlanta,  GA,  May  1999. 

5.  Optimal  Inverse  Design  using  Sensitivity  Analysis  ASME  Mechanics  and  Materials  Conference, 
Virginia  Tech,  Blacksburg,  VA,  June  1999. 

6.  Transformation  Techniques  in  Sensitivity  Computations  for  Elliptic  Systems ,  Montana  State 
University,  Bozeman,  MT,  August  1999. 
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1.  Reduced  Order  Compensator  Design  for  Distributed  Parameter  Systems ,  XX  Congresso  Na- 
cional  Matematica,  Gramado,  Brazil,  September  1997. 

2.  Sensitivity  equation  Methods  for  Optimal  Design ,  Conference  on  Multidisciplinary  Design  Op¬ 
timization,  Montreal,  Quebec,  March  1998. 

3.  A  Distributed  Parameter  Approach  to  Optimal  Sensor  Location  and  Controller  Reduction , 
DARPA  Workshop  on  Materials  Processing  and  Control,  Washington,  DC,  April  1998. 

4.  Hybrid  Sensitivity  Equation  Methods  with  Application  to  Aerodynamic  Design ,  International 
Conference  on  Distributed  Parameter  Control,  Hangzhou,  China,  June  1998. 

5.  Optimization  Based  Design  of  PDE  Systems ,  Jiao  Tong  University,  Shanghai,  China,  June 
1998. 

6.  Sensitivity  Equation  Methods  for  Optimal  Design ,  North  Carolina  State  University,  Raleigh, 
NC,  December  1998. 

7.  Computational  Methods  for  Sensitivity  Analysis  in  Optimal  Design ,  Texas  Tech  University, 
Lubbock,  TX,  February  1999. 

8.  Sensitivitivies  in  Optimal  Design ,  International  Conference  on  Optimization,  Trier,  Germany, 
March  1999. 

9.  Optimal  Control  Approaches  to  Engineering  Design:  Numerical  and  Mathematical  Issues ,  Re¬ 
gional  Meeting  of  the  AMS,  Las  Vegas,  NV,  April  1999. 


Chapter  8 

Inventions  and  Patents 


No  patents  nor  trademarks  were  applied  for  during  this  period. 
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Chapter  9 

Honors  and  Awards 


None 


